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Abstract 

In  this  paper  we  introduce  a  new  class  of  stochastic  Petri  nets  in  which  one  or  more  places  can  hold  fluid  rather 
than  discrete  tokens.  We  define  a  class  of  fluid  stochastic  Petri  nets  in  such  a  way  that  the  discrete  and  continuous 
portions  may  affect  each  other.  Following  this  definition  we  provide  equations  for  their  transient  and  steady-state 
behavior.  We  present  several  examples  showing  the  utility  of  the  construct  in  communication  network  modeling 
and  reliability  analysis,  and  discuss  important  special  cases.  We  then  discuss  numerical  methods  for  computing  the 
transient  behavior  of  such  nets.  Finally,  some  numerical  examples  are  presented. 
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1  Introduction 


One  of  the  difficulties  encountered  while  using  Petri  nets  is  that  the  reachability  graph  tends  to  be  very  large 
in  practical  problems.  Drawing  a  parallel  with  fluid  flow  approximations  in  performance  analysis  of  queueing 
systems,  we  may  define  fluid  within  a  Petri-net  to  approximate  token  movement.  Alternatively,  some  physical 
systems  have  not  previously  admitted  a  Petri  net  modeling  approach,  as  they  explicitly  contain  continuous 
fluid-like  quantities  which  are  controlled  with  discrete  logic.  This  paper  presents  a  new  methodology  for 
modeling  such  systems. 

Stochastic  fluid  flow  models  are  increasingly  used  in  the  performance  analysis  of  communications  [3, 10,  13] 
and  manufacturing  systems.  On  the  other  hand,  stochastic  Petri  nets  with  discrete  places  provide  a  useful 
framework  for  specifying  and  solving  performance  and  reliability  models  of  discrete  event  dynamic  systems 
[1,  6,  9,  17,  19].  It  is  natural  to  extend  the  stochastic  Petri  net  framework  to  Fluid  Stochastic  Petri  Nets 
(FSPNs)  by  introducing  places  with  continuous  tokens  and  arcs  with  fluid  flow  so  as  to  handle  stochastic 
fluid  flow  systems.  This  paper  extends  the  model  in  an  earlier  paper  [18]  by  allowing  the  level  of  fluid 
in  continuous  places  to  affect  the  enabling  of  timed  transitions  and  the  rates  of  fluid  flow  into  and  out 
of  continuous  places.  Rules  for  transition  enabling  and  firing  are  extended  to  reflect  the  notion  that  flow 
through  a  fluid  place  represents  token  movement. 

An  FSPN  contains  two  types  of  places:  discrete  places  containing  a  non-negative  integer  number  of 
tokens,  and  continuous  places  containing  fluid.  Transition  firings  are  determined  by  both  discrete  places 
and  continuous  places,  and  fluid  flow  is  permitted  through  the  enabled  timed  transitions  in  the  Petri  net. 
Associating  exponentially  distributed  or  zero  firing  time  with  transitions,  we  can  then  write  the  differential 
equations  for  the  underlying  stochastic  process.  We  also  provide  additional  examples  of  FSPN  usage,  and 
discuss  numerical  issues  arising  in  the  solution  of  the  underlying  dynamic  equations. 

The  main  motivation  of  this  paper  is  to  put  the  research  by  Mitra  and  his  colleagues  [3,  10,  13]  in  the 
context  of  Petri  nets,  to  make  some  extensions  and  to  investigate  the  numerical  transient  analysis  of  such 
stochastic  fluid  models. 

The  paper  is  organized  as  follows.  In  Section  2  we  develop  the  fluid  model  of  stochastic  Petri  nets  and  in 
Section  3  we  discuss  their  analysis.  Examples  and  special  cases  are  described  in  Section  4,  numerical  solution 
techniques  are  described  in  Section  5,  while  numerical  examples  are  given  in  Section  6.  The  paper  concludes 
in  Section  7. 


2  The  Stochastic  Fluid  Model 


Following  the  customary  notation  [8,  14,  16]  for  defining  Petri  nets  and  their  extensions,  we  define  a  fluid 
stochastic  Petri  net  (FSPN)  as  a  8-tuple  (J> ,T ).  V  is  a  set  of  places  partitioned  into  a 
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set  of  discrete  places  Vd  and  a  set  of  continuous  places  Vc  ■  The  number  of  fluid  places  is  F  >  0,  indexed  by 
k  =  1, 2,  •  •  • ,  F.  In  a  graphical  representation,  we  shall  depict  continuous  places  by  means  of  two  concentric 
circles.  The  set  of  transitions  7  is  partitioned  into  a  set  of  (exponentially  distributed  firing)  timed  transitions 
Te  and  a  set  of  immediate  transitions  7j.  The  set  of  directed  arcs  A  is  partitioned  into  two  subsets  Ac  and 
Ad.  Ac  is  a  subset  of  (Vc  x  TE)  U  ( TE  x  Vc)  while  Ad  is  a  subset  of  {Vd  x  T)  U  (T  x  Vd).  In  a  graphical 
representation,  arcs  in  Ac  are  drawn  as  double  lines  (to  suggest  a  pipe)  while  those  in  Ad  are  drawn  as  single 
lines. 

Let  md  —  (#p*,  i  E  Vd)  be  the  vector  of  the  number  of  tokens  in  discrete  places  and  let  x  =  (£fc,  k  £  Vc) 
be  the  vector  of  the  fluid  levels  in  continuous  places.  We  will  say  that  is  the  discrete  marking  of  the 
net.  Let  M  denote  the  number  of  discrete  markings,  which  will  be  indexed  by  the  symbols  i  and  j,  with 
i  j  —  1, 2,  •  •  • ,  M.  The  complete  state  (marking)  of  a  fluid  Petri  net  is  described  by  the  vector  m  =  (x,  md) 
where  md  is  the  marking  of  the  discrete  part  of  the  state  and  x  keeps  track  of  the  fluid  levels  in  the  continuous 
places.  Let  M  be  the  set  of  all  complete  markings  (x,  md)  and  Md  be  the  set  of  all  discrete  markings.  The 
initial  marking  is  mo  =  (xo,  md o). 

In  our  formulation  an  enabled  transition  in  7E  may  drain  fluid  out  of  its  continuous  input  places,  and 
may  pump  fluid  into  its  continuous  output  places.  The  rates  of  flow  may  be  dependent  on  the  complete 
marking  ( x ,  md).  In  a  general  formulation  of  a  stochastic  Petri  net  (embodied,  for  example,  by  SPNP  [7,  8]) 
the  conditions  for  enabling  can  be  specified  either  through  explicit  arcs  or  through  Boolean  functions  known 
as  guards.  We  will  allow  both  of  these  possibilities.  We  will  continue  to  use  the  enabling  and  firing  rules 
employed  in  SPNP  with  the  additional  possibility  of  a  guard  associated  with  a  timed  transition  being  able 
to  base  the  enabling  condition  not  only  on  the  discrete  marking  of  the  net  but  also  on  the  continuous  part. 
We  disallow  the  enabling  of  immediate  transitions  by  fluid  levels  as  this  would  lead  to  ^-dependent  vanishing 
markings,  which  cannot  be  eliminated  in  a  manner  analogous  to  that  of  GSPNs.  Thus  the  guard  function  Q 
is  defined  for  any  timed  transition  in  JE  so  that  Q  :  T  x  M  — +  {0, 1}.  For  a  timed  transition  r  £  7E,  C/(r.  rn) 
is  a  Boolean  function  that  will  be  evaluated  in  each  marking,  and  if  it  evaluates  to  irue,  the  transition 
may  be  enabled;  otherwise  r  is  disabled.  Upon  firing,  the  transition  removes  a  specified  number  of  tokens 
from  each  discrete  input  place,  and  deposits  a  specified  number  of  tokens  in  each  discrete  output  place. 
The  basic  extension  we  have  made  from  [18]  here  is  that  fluid  levels  in  continuous  places  can  change  the 
enabling/disabling  of  timed  transitions  in  the  discrete  part  of  the  net.  A  discrete  marking  md  is  said  to  be  a 
vanishing  marking  if  one  or  more  immediate  transitions  are  enabled  by  it;  otherwise  it  is  a  tangible  marking. 

The  firing  rate  function  I  is  defined  for  timed  transitions  JE  so  that  T  \  7  E  x  M  —>  IR+.  Thus  if  a 
timed  transition  r  is  enabled  in  (tangible)  marking  m,  it  fires  with  rate  7(r,  m).  Note  once  again  that  in 
[18],  these  rates  were  not  allowed  to  depend  upon  the  fluid  levels  but  now  we  do  allow  the  firing  rates  to  be 
dependent  on  fluid  levels. 

As  in  [18],  the  weight  function  W  is  defined  for  immediate  transitions  7}  so  that  W  :  7 j  x  Md  — *  IR+. 
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Thus  if  an  immediate  transition  i  is  enabled  in  (a  vanishing)  marking  md,  it  fires  with  probability 

_ W(i,  rnd) _ 

53  W(0,md) 

OeTT  enabled  in  m<* 


Next  we  describe  the  evolution  of  the  continuous  part  of  the  marking.  The  flow  rate  function  1Z  is  defined 
for  the  arcs  connecting  a  continuous  place  and  a  timed  transition  so  that  ft  :  Ac  x  M  ^  IR+  U  {0}.  Thus 
when  the  FSPN  marking  is  mt  €  M  at  time  t ,  fluid  can  leave  place  k  £  Vc  along  the  arc  ( k ,  r)  £  Ac  at  rate 
H((i and  can  enter  the  continuous  place  k  at  rate  k),mt)  along  the  arc  (r,  k)  £  Ac  for  each 
r  £  JE  that  is  enabled  in  mt.  The  instantaneous  rate  at  which  fluid  builds  in  a  place  k  £  Vc  at  time  t,  in 
marking  mt ,  is  then  given  by 

r*(m*)=  53  ft((r,fc),mt)-  53  n{(k,r),mt). 

t£Te  enabled  in  mt  r^TE  enabled  in  mt 

We  require  that  for  every  discrete  marking  md  and  arc  (r,  k),  the  rate  k),  (x,  md))  be  a  “nice”  function 
of  x,  e.g.,  it  is  piecewise  continuous.  Observe  that  since  mt  contains  continuous  levels  x,  rk(mt )  may  change 
as  a  function  of  t  even  if  the  discrete  part  of  mt  does  not  change.  Once  again  we  have  extended  the  definition 
in  [18]  by  allowing  these  rates  to  be  dependent  on  the  fluid  levels. 

Now  let  Xfc(t)  be  the  fluid  level  at  time  t  in  a  continuous  place  k  e  Vc.  We  assume  that  there  is  an  upper 
bound  on  the  fluid  content,  that  is,  Xk (t)  <  Bk  for  all  t  >  0.  If  there  is  no  such  upper  bound,  we  set  Bk  to 
oo.  Then  the  sample  path  of  Xk(t)  satisfies  the  differential  equation 

IMro<)]+  if 

rk(mt)  if 

0  if 

In  the  case  Xk(t)  —  0  and  rk(mt)  <  0,  we  set  the  actual  rate  equal  to  zero  (denoted  by  [rfc(mt)]+  = 
max(rfc(raf),  0))  in  order  to  maintain  Xk(t)  >  0.  In  the  case  that  Xk(t)  =  Bk  and  rk(mt)  >  0,  we  set  the 
actual  rate  equal  to  zero  (denoted  by  [rk(mt)]~  =  min( rk(mt),  0))  in  order  to  maintain  Xk(t)  <  Bk.  For  the 
explantion  of  the  remaining  cases,  we  refer  the  reader  to  [10],  Section  II.  The  key  observation  (for  the  fourth 
case)  is  that  a  sign  change  from  +  to  -  in  rk(mt)  at  mt  will  “trap”  Xk(t)  to  be  constant.  Finally,  let  Md(t) 
be  the  discrete  marking  at  time  t. 

In  the  next  section  we  study  the  joint  process  (X(t),  Md(t)),  where  X(t)  =  [Xfc(t),  k  €  Vc]- 


Xk(t)  =  0 
Xt(t)  =  Bk 

0  <  Xk(t)  <  Bk  and  rk(mt-.)rk(mt+)  >  0 
0  <  Xk(t)  <  Bk  and  rk(mt-)rk(mt+)  <  0 
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3  Analysis 


Recall  that  Xk(t)  is  the  fluid  level  in  the  kth  continuous  place  at  timet.  The  reachability  graph  corresponding 
to  the  discrete  part  of  the  net  gives  rise  to  a  stochastic  process  that  is  a  Markov  process  with  the  state  space 
M.  In  the  special  case  that  the  discrete  part  of  the  net  is  not  affected  by  the  fluid  levels,  discrete  part  of 
the  net  gives  rise  to  a  CTMC[1].  Let  5  be  the  discrete  state  space  and  let  Q(x)  =  [qij  (£)]  be  the  matrix 
of  transition  rates  derived  from  the  firing  rate  function  T  of  section  2.  S  corresponds  to  the  set  of  tangible 
discrete  markings  in  Md •  For  every  x  and  place  k  E  Vc  define  the  diagonal  matrix  Rfc(x)  =  diag(rk(x,  md)), 
m,d  €  5). 

Define  the  distribution  function  H(t,x,m.d)  =  P  (X(t)  <  x,  and  Md{t)  =  md)  and  let  H(t,x )  = 
[H(t,  x,  md),  md  £  5]  be  a  row  vector. 


To  begin  with,  we  assume  that  the  rate  functions  Rk(x)  are  differentiable  functions  of  x  and  transition 
rates  Q(F)  are  piecewise  right  continuous  functions  of  x.  We  also  assume  that  the  capacities  of  the  continuous 
places  are  infinite.  Under  these  assumptions,  it  can  be  shown  that  (X(t),  Md(t))  has  a  density  h(t,x,md) 
for  all  x  with  non-negative  components  and  md  €  S  and  has  a  probability  mass  c(t,x,  md)  if  x  has  at  least 
one  component  equal  to  0,  where 


c(t,  x,  md) 


P(Xk(t)  =  0  if  Xk  =  0,  Xkjt)  €  (xk,  Xk  +  dxk)  if  Xk  >  0  VjQ 

rifc:Xfc>0  dXk 


Let  h(t ,  x)  denote  the  corresponding  row  vector  of  h(t,  x,  md). 

The  next  theorem  gives  the  coupled  system  of  partial  differential  equations  satisfied  by  h  and  c.  These 
equations  describe  the  transient  behavior  of  the  FSPN. 


Theorem  1: 

The  equations  satified  by  h  are 


dh  v  ■> 
dt  dxk 

k 

c(t,x,md) 


— c(t,x,md)+  ^2  h(t,x,md)rk(x,md)+ 

™  k:xk=0 

^2 


k:xk>0 


=  hQ(x) 

=  0 

if,  for  any  k,  Xk  =  0  and  rk(x,md)  >  0 

=  ^  ^ c(t,  X,  i)qi,mjiP) 

ies 

if  Xk  =  0  =>  rk(x,md)<  0  Vfc 


(2) 

(3) 


(4) 
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Proof: 


A  rigorous  proof  can  be  obtained  by  using  the  same  technique  as  in  [12].  Here  we  provide  an  intuitive 
proof  for  the  ft-equation. 

Assume  for  simplicity  only  one  fluid  place  and  ri(*,i)  >  0.  Consider  a  portion  of  the  function  h(t,x,i) 
at  one  instant  in  time  in  the  vicinity  of  the  location  x  =  x\.  We  consider  a  cell  surrounding  this  point  whose 
left  and  right  boundaries  are  located  at  x _  and  x+  respectively.  Let  the  discrete  values  ‘h(t,xi,i)  and  qij(xj) 
represent  the  mean  values  of  h  and  qij  respectively  within  the  cell.  This  situation  is  illustrated  in  Figure  1. 


h(t,x,i) 


h(t,x,j ) 


x 


xl 


x_ 


♦ 


xl 


I _ 

Ax 


J 


Figure  1:  Derivation  of  FSPN  equations 

Probability  mass  must  be  conserved.  We  may  therefore  derive  a  balance  equation  for  the  change  in 
probability  mass  inside  the  cell: 


change  in  probability 

total  mass 

total  mass 

mass  in  cell 

entering  cell 

leaving  cell 

Probability  mass  enters  the  cell  at  the  location  x _  at  the  rate  r(x-,i)  and  leaves  at  at  the  rate  r(x+1i). 
In  addition  the  cell  gains  probability  from  the  corresponding  cell  of  each  other  discrete  marking  mj ,  j  ^  i 
at  the  rate  qji(xi)  and  loses  probability  to  them  at  a  total  rate  of  -qn(x\ ).  The  balance  equation  for  the  cell 
probability  mass  becomes 

A -  h(t,  x-,i)n(x- ,  i )  -  h(t,  x+,i)r1(x+,i)  +  Ax(]T qji(xi)h(t, xhj)  +  qu(xi)h(t,  x{, «')).  (5) 

The  equation  is  identical  for  the  case  ri(x ,  i)  <  0.  Dividing  the  vector  form  of  equation  (5)  by  Ax  and  taking 
the  limit  as  Ax  — ►  0  we  obtain 

g+2*|M  =  SQ(f). 

The  argument  generalizes  easily  to  the  case  of  more  than  one  fluid  place,  yielding  equation  (2).  The  cell  is 
in  general  a  hypercube  of  sidelength  Ax. 
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The  boundary  conditions  follow  from  the  fact  that  h(t,x,  i )  is  a  probability  density  function. 


Theorem  2: 

The  equations  for  the  cumulative  probability  distributions  H  are 

f + £  = B «<*>  -  f  ■■■['  dtl -*CF 


with  the  boundary  conditions 


H(t,x,i)  =  0  at  xk  =  0,  if  rk(x,i)  >  0 


r  dH  n 
lim  « —  =  0. 

oo 


Proof: 


Using  the  abbreviations 


. . .  dx f 


and  noting  that 


px  rXF  rx  i 

dx  =  /  •  •  •  /  . . .  dxf 

Jo  JO  Jo 

H(t,x,i)  =  c(t,x,i)+  [  h(t,  x,  i)dx 
Jo 


:  substitute  into  Equation  (2)  and  integrate  with  respect  to  x: 


which  yields  (6).  Equation  (4)  is  implicitly  contained  in  Equation  (6),  along  with  the  boundary  condition 

(7). 

The  boundary  condition  (7)  follows  from  the  observation  that  the  fluid  level  in  place  k  cannot  remain 
at  0  for  a  positive  amount  of  time  in  state  if  the  net  rate  ruj)  >  0.  The  boundary  condition  (8) 
comes  from  the  observation  that  H  represents  a  probability  distribution  with  respect  to  each  xk,  and  must 
therefore  approach  an  asymptotic  value  xk  tends  to  infinity  (or  reaches  its  maximal  value). 


This  method  of  deriving  partial  differential  equations  that  represent  conservation  laws  is  well  established 
in  computational  fluid  dynamics,  where  it  is  known  as  the  “Finite  Volume”  technique  [15]. 
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The  assumptions  on  Rfc(x)  can  be  relaxed  and  we  can  allow  Kk(x)  to  be  piecewise  differentiable  and 
fluid  places  to  have  finite  capacities.  This  may  introduce  non-zero  probability  masses  at  the  inter-region 
boundaries  and  will  need  to  be  explicitly  accounted  for. 

The  domain  of  Equations  (2)  and  (6)  is 


0  <  xk  <Bk 

0  <  t  <  oo, 

although  in  practice  we  will  only  be  interested  in  the  finite  domain  0  <  t  <  tmax ,  0  <  Xk  <  min {xmax,  Bk} 
(where  xmax  is  finite  when  Bk  =  oo). 

The  initial  conditions  for  Equation  (6)  are 

H( 0,  x ,  i)  =  1  if  i  =  rrido  and  x  >  x0 

=  0  otherwise. 

The  initial  conditions  for  Equation  (2)  are 

c(0,£,i)  =  <5(m0) 


where  S  is  the  delta  function. 


Now  suppose  the  following  limits  exist 

F(x)  =  lim  H(t,x). 

t — ►  oo 

Then  from  Theorem  2,  we  see  that  the  steady-state  distribution  Fix)  obeys  the  following  system  of  differential 
equations 

£  £;  (**>**<*>) = Ji~  <9> 
k 

with  the  normalization  condition  lim  F(x)e  =  1,  where  e  is  a  column  vector  of  all  l’s. 

X — *■  OO 

We  note  that  the  steady-state  distribution  F  exists  when 

lim  Ti(x)rk(x,i)  <  0  =  l  ...F 

art— » oo 

i 

where  7r(x)  is  the  solution  of 


7r(f’)Q(x)  =  0 

53  *■*'(*)  = 1 
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3.1  FSPN  with  a  Single  Continuous  Place 

In  the  special  case  of  a  single  continuous  place,  Equation  (9)  reduces  to: 


±(F(x)R(x))  =  F(x)Q(x)  -  £  F^j£  dx. 

In  the  following  subsections,  we  consider  three  special  cases  of  the  Equation  (10). 


(10) 


3.1.1  Constant  Case 

In  the  special  case  that  R(ar)  and  Q(s)  are  both  independent  of  x,  following  [3],  solution  of  such  an  equation 
is  of  the  form  F(x)  =  heXx  where  h  is  a  row  vector  and  A  is  a  scalar.  Substituting  in  (10)  we  have 


ft(AR-Q)  =  0. 


(11) 


If  a  non-zero  h  is  to  satisfy  the  above  equation,  we  must  have  det(AR  —  Q)  =  0.  The  number  of  solutions 
of  det(AR  -  Q)  =  0  equals  the  number  of  non-zero  diagonal  elements  of  R(x) .  Let  these  solutions  be  denoted 
by  Ai,  A2  •  •  ■ ,  A*,.  Let  hi  be  the  solution  to  A,-( A*R  -  Q)  =  0.  Then  the  general  solution  to  (10)  is  given  by 

k 

F{x)  =  J2a^eXiX  (12) 

2  =  1 

where  the  scalars  a;  need  to  be  determined  from  the  boundary  conditions  and  the  boundedness  of  F(x).  It 
is  known  that  the  number  of  A,’s  with  positive  real  part  equals  the  number  of  negative  diagonal  entries  of 
R.  The  coefficient  a,  corresponding  to  an  eigenvalue  A;  with  Re(\t)  >  0  must  be  zero  in  order  to  maintain 
boundedness  of  F(x).  The  remaining  coefficients  a,  are  uniquely  determined  by  the  boundary  condition 
F(0,  m)  =  0  if  r(m)  >  0. 

Now  if  R(x)  and  Q(x)  are  both  piecewise  constant  functions  of  x,  we  can  apply  the  above  procedure 
for  each  different  segment  and  piece  the  individual  solutions  together  [10].  In  the  general  case,  we  can  use 
numerical  solution  methods  for  linear  odes  that  are  available.  We  refer  here  to  explicit  methods  such  as 
RKF-45  or  implicit  methods  such  as  implicit  Runge-Kutta  [4]  or  TR-BDF2  [5]  and  so  on. 


4  Examples 

Next  we  examine  a  number  of  examples  to  illustrate  the  modeling  power  of  FSPNs. 
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Figure  2:  FSPN  Model  of  Statistically  Multiplexed  Network  Switch 


4.1  Statistical  Multiplexing  of  a  Network  Switch 

For  our  first  example,  we  recast  the  problem  studied  in  [10]  as  a  FSPN.  We  have  K  sources  of  ATM  cells. 
Each  source  emits  both  high  priority  and  low  priority  cells;  the  arrival  rate  of  each  depends  on  the  state  of 
a  two-state  Markov  chain,  e.g.,  if  the  source’s  Markov  chain  is  in  state  $  E  [1,2],  then  the  rate  at  which  high 
priority  cells  is  produced  is  A^,  and  the  rate  at  which  low  priority  cells  is  produced  is  A^.  All  cells  are 
delivered  to  a  single  switch,  with  buffer  capacity  B.  All  low  priority  cells  are  discarded  when  the  buffer  level 
is  greater  than  some  B%  <  B. 

The  FSPN  for  this  problem  is  illustrated  in  Figure  2.  The  number  of  tokens  in  place  pi  reflects  the 
number  of  sources  in  state  i,  i  =  1,2.  The  rate  at  which  sources  change  from  state  1  to  2  (alt.,  2  to  1) 
is  #piAi  (alt.,  #P2^2  ).  Letting  rfigh  and  r\ow  denote  the  rate  at  which  high  and  low  priority  cells  are 
generated  by  a  source  in  state  i,  the  aggregate  arrival  rate  of  fluid  to  the  buffer  from  sources  in  state  i  is  a 
function  of  the  discrete  marking  : 

fi(#Pi,#P2)  =  #Pi(r?9h+rl™). 

The  overall  arrival  rate  as  a  function  of  the  discrete  marking  is 

7(#Pl  >  #P2)  =  /l(#Pl ,  #P2)  +  /2(#Pl ,  #P2). 

The  fluid  level  in  the  place  loss  reflects  the  total  cell  loss  since  time  0  (like  the  cell  buffer,  it  is  initially 
empty).  The  rate  c(x)  at  which  cells  are  successfully  switched  out  of  the  buffer  with  level  x  is  c(x)  =  c  if 
x  >  0,  and  c(x)  =  0  if  x  =  0.  Of  more  interest  is  the  function  T(#Pi,  #P2, x)  describing  the  rate  of  cell  loss 
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Figure  3:  FSPN  for  the  Machine  Breakdown  Model 


(both  into  and  out  of  the  transition  so  labeled): 


L(#Pi,#P2>x) 


'  [7(#Pi.#P2)-c]+  ifO<x<Bt 

#Pir[™  +  #PA0W  +  [#Pir\i9h  +  #p2r \%sh  -  c]+  if  Bt  <  x  <  B 

#Pir'r  +  #P2 r‘20W  +  [#Pirf 5/1  +  #V2rh7igh  -  c]+- 

[c-#p1r^ft+#p2r2ft^]+  iix  =  Bt 


When  x  <  Bt,  low  priority  cells  are  not  discarded  automatically,  and  any  loss  is  the  difference  between  the 
aggregate  arrival  rate  and  the  service  rate.  Low  priority  cells  are  dropped  automatically  when  Bt  <  x  <  B, 
and  further  loss  may  be  due  to  an  inability  of  the  server  to  keep  up  with  the  aggregate  high  priority  arrival 

rate. 


4.2  Machine  Breakdown  Model 

In  the  previous  example,  fluid  levels  have  no  effect  on  the  behavior  of  the  discrete  portion  of  the  FSPN.  The 
next  example  shows  how  fluid  levels  may  affect  the  firing  rate  of  discrete  transitions.  Consider  a  system  with 
N  statistically  identical  and  independent  components,  each  with  a  failure  rate  that  depends  on  the  overall 
system  load  [11];  the  repair  rate  is  p.  Work  arrives  to  the  system  at  a  constant  rate  r  and  is  completed  at 
rate  d  per  functioning  machine.  Work  here  is  considered  to  be  a  non-negative  real  quantity.  We  model  this 
system  as  an  FSPN  shown  in  Figure  3. 

The  model  has  two  discrete  places  (p„  and  pd),  one  continuous  place  and  three  timed  transitions  (tf,  tr 
and  t always)-  The  number  of  tokens  in  pu  models  the  number  of  functioning  machines  (with  the  initial  value 
N)  while  the  number  of  tokens  in  pd  is  the  number  of  failed  machines  undergoing  repair.  The  firing  of  the 
transition  tf  represents  a  machine  failure;  given  positive  failure  rates  Aj  and  A2  and  positive  scaling  factor 
a,  the  load-dependent  firing  rate  of  this  transition  with  fluid  level  x  is  taken  to  be 

f(#Pu,x)  =  #pu(Ai  +  A2(1.0  -  e-a*)). 
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Off 


Figure  4:  FSPN  model  of  an  alternating  switch 


Repair  of  a  machine  is  modeled  by  the  firing  of  tr]  its  firing  rate  is  p  times  the  number  of  tokens  in  pd -  The 
transition  taiways  is  always  enabled  and  it  continuously  pumps  fluid  at  rate  r  into  the  continuous  place.  The 
rate  of  firing  of  this  transition  can  be  chosen  to  be  any  positive  value.  Whenever  transition  tj  is  enabled,  it 
drains  fluid  from  the  continuous  place  at  rate  d  times  the  number  of  tokens  in  place  pu . 


4-3  Alternating  Switch 

In  the  next  example  we  model  a  communcation  switch  with  two  arrival  streams.  One  stream  is  bursty,  with 
a  high  arrival  rate  rhigh  when  active.  It  is  described  as  an  alternating  process,  that  is  on  for  an  exponential 
period  of  time  with  rate  Xon ,  and  off  for  an  exponential  period  of  time  with  rate  X0j j .  The  other  stream  is 
slow,  but  constant  (riow).  The  switch  services  workload  from  either  stream  at  rate  p. 

The  switch  is  designed  to  allocate  exponentially  distributed  time-slots  to  each  stream,  with  rates  A  high 
and  Xiow .  The  time  slots  alternate — the  fast  stream  is  given  a  slot,  then  the  slow  stream,  and  so  on.  The 
FSPN  modeling  this  switch  is  illustrated  in  Figure  4. 
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5  Numerical  Transient  Analysis 


5.1  Discretization 


We  choose  to  solve  the  equations  for  the  probability  distribution  (6)  rather  than  the  density  equations 
(2),  since  the  latter  would  require  the  numerical  treatment  of  a  delta  function  used  to  describe  the  initial 
condition.  We  describe  here  the  case  R(x)  =  const,  Q(x)  =  const. 


We  perform  a  semi-discretization  of  (6)  in  the  coordinate  directions  xk,  using  a  first-order  upwind  method. 
In  the  nomenclature  of  section  3,  the  upwind  scheme  approximates  values  at  the  boundaries  of  the  cell  by 
those  at  the  neighbouring  grid  points  in  the  upstream  direction,  relative  to  the  flow  direction  r.  We  choose  a 
finite  domain  0  <  xk  <  xmax  and  use  G  equidistant  grid  points.  The  grid  spacing  (mesh  size)  Ax  is  therefore 
equal  to  xmax/(G  -  1).  Note  that  in  order  that  the  boundary  condition  (8)  be  approximately  satisfied, 
xmax  must  be  chosen  to  be  sufficiently  large.  We  use  the  notation  rk,iu...,iF,i  =  rk(hAx, . . lFAx,  i), 
=  qij(hAx,...,lFAx)  and  HiM . ,F(t )  =  H(t,hAx, . . . , lFAx,  i)  .  The  upwind  discretization  is 

given  by 


— — H(t , . . . ,  Xk  = 
OXk 


■t  ■■■(*)  -  </-D»....(0)  if  . . .  £  0 

&  (Hi (i+1  )*,...(*)  -  Hi 7fcl  .(0)  if  rk,...,ih,...,i  <  0 


(13) 


We  obtain  from  the  semi-discretization  the  linear  system  of  ordinary  differential  equations 

^-  =  HQ  (14) 

dt 

where  the  unknown  vector 

H  =  (jfiTi.o, . 0F>  f.  •••>  Hm,G-U . <?— If) 

is  obtained  by  a  lexicographic  ordering  of  the  unknowns  by  fluid  place  and  by  the  discrete  marking. 

The  matrix  Q  is  given  by 

Q  =  Q  +  W 

where  Q  represents  the  CTMC  of  the  discrete  part  of  the  net,  and  W  the  discretization  of  the  space 
derivatives  multiplied  by  the  flow  rates. 


The  matrix  Q  is  given  by 


Q  n  •••  Qim 

q  =  ;  .  ; 

Qmi  •  •  •  Qmm 
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where 


Qij  — 


#17,01,.,., Of 


Note  that  the  matrix  Q  is  a  CTMC  matrix.  The  block-diagonal  matrix  W  has  the  form 


W  =  - 


_1_ 

Ax 


Wi 


0 


0 

WM 


where  Wi  €  IR<?FxgF  represents  the  discretization  of  the  second  term  of  equation  (6)  for  the  i-th  discrete 
marking  using  the  upwind  method  (13).  Each  Wi  is  a  sparse,  banded  matrix  in  which  each  column  is  either 
zero  or  has  a  positive  main  diagonal  coefficient  and  non-positive  off-diagonal  coefficients  containing  the  rate 
values  ....  In  addition,  the  main  diagonal  coefficient  is  the  negative  sum  of  the  other  entries  in  its  column. 
Each  Wi ,  and  thus  also  W,  therefore  represents  the  transpose  of  a  CTMC. 


The  linear  system  of  ordinary  differential  equations  (14)  thus  defined  may  be  integrated  by  any  standard 
method.  In  the  numerical  examples  in  section  6  we  will  use  the  Forward  Euler  scheme 

H(t  +  At)&H(t)  +  AtH(t)Q. 


Note  that  for  this  scheme  the  discretization  mesh  sizes  must  satisfy 


max  rjg  i  At  <  Ax 

fc,?!,...,  lp,i  ’  ’ 


(15) 


in  order  that  the  integration  be  stable. 


5.2  Numerical  Integration  in  a  Transformed  Domain 

In  order  to  avoid  the  difficulty  of  solving  equation  (6)  in  an  arbitrarily  large  and  a  priori  unknown  domain, 
we  can  perform  the  coordinate  transformation 

Vk  =  1  -  (16) 

which  maps  the  infinite  interval  Xk  —  [0, 00]  to  the  finite  interval  yk  =  [0, 1].  Equation  (6)  then  becomes 

f + ?  £r‘w<i  -  »> = *«*>  -  j:  Bwdy  (n) 

for  H(i,y).  The  boundary  and  initial  conditions  remain  unchanged. 

This  minor  modification  to  the  equations  makes  their  solution  significantly  easier,  since  the  decision 
where  to  place  xmax  must  no  longer  be  made,  and  the  danger  of  an  inappropriate  choice  is  avoided.  The 
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coordinate  transformation  has  the  additional  advantage  of  compressing  the  infinite  interval  of  the  large  x 
values,  where  in  many  cases  the  solution  shows  virtually  no  structure,  into  a  small  space.  Furthermore, 
we  can  obtain  numerical  values  at  x  =  oo,  which  correspond  to  the  probabilities  for  the  discrete  markings, 
whereas  in  the  untransformed  case  these  must  be  approximated  by  values  obtained  at  x  =  xmax. 

5.3  Problem  Size 

Recall  that  G  denotes  the  number  of  gridpoints  in  each  dimension  Xi,  F  the  number  of  fluid  places,  and  M 
the  number  of  discrete  markings.  The  number  of  time-steps  to  be  integrated  is  T. 

The  computational  complexity  for  the  solution  is  0{TMGF )  floating  point  operations,  since  for  each  of 
T  timesteps  we  must  increment  each  solution  value  in  an  F-dimensional  grid  of  sidelength  G  for  each  of 
M  discrete  markings  using  0(1)  operations.  Note  that  for  an  explicit  integration  method  such  as  Forward 
Euler,  because  the  condition  (15)  must  be  satisfied,  an  increase  in  G  must  ultimately  be  accompanied  by  a 
proportionate  increase  in  T . 

The  storage  requirements  of  the  algorithm  are  at  least  4 MG^  bytes  since  for  each  of  M  discrete  markings 
we  must  store  an  F- dimensional  grid  of  floating  point  numbers  with  sidelength  G.  Solutions  at  successive 
time-steps  can  be  overwritten. 

Since  the  sidelength  of  the  grids  G  may  typically  be  of  the  order  of  50  or  more  when  a  simple  discretization 
is  used,  we  see  that  this  would  seriously  limit  the  size  of  the  FSPNs  that  can  be  solved.  Future  work  must 
therefore  include  strategies  for  reducing  the  amount  of  memory  needed  to  represent  the  function  H. 

5.4  Time  Integration  by  Randomization 

Randomization  is  a  numerical  method  widely  used  for  the  solution  of  systems  of  ODEs  of  the  form  (14).  It 
has  the  advantages  of  high  numerical  stability  and  low  roundoff  error,  a  priori  specification  of  absolute  error 
tolerance  requirements,  and  is  in  addition  often  found  to  be  faster  than  numerical  integration  schemes.  The 
method  has  superior  roundoff  error  behavior  when  the  matrix  P  =  (/  +  yQ)  has  all  non-negative  entries, 
where  A  >  max|Q*i|.  This  is,  for  example,  the  case  for  a  CTMC. 

As  the  following  Theorem  shows,  the  semi-discretized  FSPN  equation  (14)  also  has  this  property,  indi¬ 
cating  that  randomization  may  be  the  method  of  choice  for  computing  transient  solutions. 

Theorem  3: 

For  the  matrix  P  =  (1  +  AQ)  where  Q  is  the  matrix  of  equation  (14)  holds 

P ij  >  0  Vi,  i 
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Figure  5:  Numerical  Solution  of  Breakdown  Model 


Proof: 

It  is  sufficient  to  show 

Q«j  W«  <  0,  and  Qfj,  Wy  >0  i  /  j  (18) 

“Qn  >  Q ij,  Q ji  i  #  j  (19) 

We  have 

Qn  ?  Wii  <  0,  Qz; ,  W ij  >0  i  ±  j 

Qu  =  -X>y,  Wii  =  -^Wi7- 

since  Q  is  a  CTMC  and  W  represents  an  upwind  discretization.  Equations  (18)  and  (19)  follow  directly 
from  Q  =  Q  +  W. 


6  Numerical  Examples 

6.1  Machine  Breakdown  Model 

For  our  first  illustration  of  the  behaviour  of  an  FSPN  we  choose  the  model  of  section  4.2.  We  consider  the 
case  of  one  processor  only,  and  parameter  values  ofAi  =  2,  A2  =  0,  //  =  3,  r  =  1  and  d  =  2.  In  this  case, 
both  R  and  Q  are  independent  of  x .  We  solve  the  equations  for  the  distribution  (6)  in  the  range  0  <  t  <  4, 
0  <  x  <  4.  We  discretized  with  stepsizes  Ax  =  1/64  and  Atf  —  l.Oe  -  4.  The  numerical  results  for  this 
problem  are  shown  in  Figure  5. 

For  this  simple  case,  it  can  be  shown  via  Laplace  transforms  that  the  steady  state  solution  is  given  by 

H(x,UP)  =  0.4(1  -  e~x) 
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H(DOWN) 


Figure  6:  Numerical  Solution  of  Modified  Breakdown  Model 


H(x,DOWN)  =  0.2(3  -  2e~x) 

The  transient  solution  values  obtained  at  t  =  4  agree  closely  with  these  results. 

Now  we  modify  the  example  to  allow  dependency  of  the  matrices  R  and  Q  on  the  fluid  level  X.  First  we 
extend  the  model  to  contain  a  backup  server  which  is  only  used  when  the  primary  server  is  down  and  the 
load  level  reaches  a  value  of  0.5.  We  thus  have 

{2  md  =  UP 

0  mi  —  DOWN,  x  <  0.5 
2  md  =  DOWN,  x  >  0.5 

The  dependency  of  Q  on  x  is  obtained  by  setting 

A  =  Ai  +  A2(l  —  e~ax) 

choosing  a  =  1,  Ai  =  1  and  A2  =  2.  The  results  of  this  computation  are  depicted  in  Figure  6.  Note  the 
discontinuity  and  change  of  slope  at  x  —  0.5,  when  the  backup  server  is  started. 

Figure  7  shows  the  results  for  the  unmodified  breakdown  model,  using  the  equation  in  transformed 
coordinates  (17).  Here  integration  until  t  -  24  has  been  performed.  Note  that  the  solution  obtained  for 
H(DOWN)  at  t  =  24  appears  as  a  straight  line  through  the  origin  with  gradient  0.4,  corresponding  to  the 
analytic  solution  of  Equation  (20). 


6.2  ATM  Switch  Model 


We  now  consider  the  multiplexed  network  switch  model  considered  in  section  4.1.  We  set  the  number  of 
sources  I<  to  one,  parameters  Bl  =  3000  and  B  =  6000.  High  priority  packets  arrive  at  the  rate  6  x  103  and 
low  priority  packets  at  a  rate  of  4  x  103  per  second  when  the  source  is  in  state  1,  and  at  rates  4  x  103  and 
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H(UP)  H(DOWN) 


Figure  7:  Probability  distributions  for  breakdown  model  in  transformed  coordinates.  Left:  Server  up;  Right: 
Server  down. 


2  x  103  per  second  respectively  when  the  source  is  in  state  2.  The  switch  can  process  both  types  of  packet 
at  a  rate  of  5  x  103  per  second.  The  exponentially  distributed  rates  of  change  between  high  and  low  priority 
packet  generation  are  A*  =0.5  and  A2  =  0.5.  The  results  for  the  number  of  packets  in  the  buffer  are  shown 
in  Figure  8.  The  left  and  right  results  at  t  =  10  qualitatively  match  those  of  Elwalid  and  Mitra  [10],  Figure 
2,  right  and  center,  respectively.  An  appropriate  choice  of  parameters  also  yields  a  result,  not  illustrated 
here,  which  is  similar  to  [10],  Figure  2,  left. 


7  Conclusion 


We  have  defined  a  new  class  of  stochastic  Petri  nets  by  introducing  places  with  continuous  tokens  and  arcs 
with  fluid  flows.  This  new  class  of  fluid  stochastic  Petri  nets  (FSPNs)  should  be  useful  in  modeling  stochastic 
fluid  flow  systems,  and  may  also  be  useful  in  modeling  processes  that  control  physical  systems.  Our  model 
formulation  permits  the  discrete  and  continuous  parts  to  affect  each  other,  endowing  FSPNs  with  the  ability 
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to  both  control  the  fluid  flow,  and  have  the  discrete  control  decisions  be  affected  by  observed  fluid  flow. 
We  have  provided  formal  definition  of  FSPNs  and  developed  the  rules  for  their  dynamic  evolution.  We 
have  derived  coupled  systems  of  partial  differential  equations  for  the  transient  and  the  steady-state  behavior 
of  FSPNs.  Spectral  representation  of  the  FSPNs  with  a  single  continuous  place  can  be  adapted  from  the 
literature  on  stochastic  fluid  flow  models.  We  have  presented  a  number  of  examples  illustrating  the  modeling 
power  of  FSPNs,  have  considered  issues  arising  in  the  numerical  solution  of  the  dynamical  equations,  and 
have  provided  numerically  solved  examples. 
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